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ABSTRACT 

According to the linearized solutions for thermal 

blooming, the density perturbations become infinite (i.e. 

"catastrophic" defocusing) as the Mach number approaches 

unity. However, the nonlinearities in the transonic equations 

cutoff the trend to infinity, and the values of the flow 

perturbation quantities are finite. The nonlinear equations 

with heat addition are transformed into simple linear 

algebraic equations through the specification of the 

streamline geometry in the heat release region. At a Mach 

number of unity, streamtube area variation was found to be 

directly proportional to the change in total temperature. 

A steady, two-dimensional mixed flow solution has been found 

for the transonic thermal blooming problem. The solution 

for the density perturbations within a laser beam at a Mach 

number of precisely unity is given. For a Gaussian beam 

7 2 

with an intensity of 3.333x10 Watts/m and an atmospheric 
absorption of 8.0x10 ^ cm ^ the maximum fractional density 
perturbation is 1.028x10 . The transonic thermal blooming 

problem does not pose as serious a problem as previously 
anticipated. 
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NOTATION 



English 

A area, square centimeters 

a speed of sound, meters/second 

b function, dimensionless 

C line contour around integration region R 

c speed of light, meters/second 

Cp specific heat at constant pressure. 

Joules/kilogram mass - degree Kelvin 

Cv specific heat at constant volxime. 

Joules /kilogram mass - degree Kelvin 

f function, dimensionless 

G Green's function, dimensionless; also specified 

function, dimensionless 

H Heat quantity, dimensionless 

h specific enthalpy. Joules/kilogram mass 

h(x,y) heat addition function. Watts/cubic centimeters 

I laser beam intensity. Watts/square centimeter; 

also integral function representation 

I(x) unit function; zero for x < 0 and unity for x ^ 0 

k Gladstone-Dale constant, cubic meters/kilogram mass 

L scaling factor, dimensionless; length, meters; 

and functional representation of normalized velocity, 
dimensionless 

£ characteristic length (subscript 1 for x-direction 

and 2 for y-direction) , centimeters 

M Mach number, dimensionless 

N Number of waves, dimensionless 

n normal coordinate, centimeters; also index of 

refraction, dimensionless 
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p static pressure. Newtons/square meter 

Pr Prandtl number, dimensionless 

Q heat source, Watts/cvibic centimeter 

q heat source. Watts/square centimeter; also, 

strength of line heat source. Watts/centimeter 

R universal gas constant. Joules/kilogram mass - 

degrees Kelvin ; radius of curvature, centimeters; 
and region of integration 

r distance, centimeters 

Re Reynolds number, dimensionless 

S line contour around shock waves; also, distance 

along a characteristic, centimeters 

s specific entropy. Joules/kilogram mass - degree 

Kelvin ; also streamline coordinate 

T static temperature, degrees Kelvin 

t maximum streamtube thickness, centimeters 

U velocity vector, meters/second 

U,u speed (x-direction) , meters/second 

V speed (y-direction) , meters/second 

w laser beam spot size, centimeters; 

also speed, meters/second 

X coordinate, centimeters 

y coordinate, centimeters 

Z streamtube shape, centimeters 

Greek 

a atmospheric absorption, 1/centimeters 

2 1/2 

3 variable, dimensionless; (1-M ) for subsonic 

flow, dimensionless; 

2 1/2 

and (M - 1) ' for supersonic flow, dimensionless 

Y ratio of specific heats, dimensionless 
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A small change, dimensionless 

6 slope, dimensionless; also delta (impulse) function 

e dimensionless small quantity, 0 < e << 1 

n absolute viscosity, kilogram mass/meter-second; 

also coordinate (y-direction) , centimeters 

6 flow angle, radians 

K thermal conductivity. Watts/centimeter - 

degree Kelvin 

X wave length, cent ime t er s 

y Mach angle defined by arcsin(l/M) 

V inverse radius of curvature, 1/centimeter 

C coordinate (x-direction) , centimeters; also 

transonic similarity parameter, dimensionless 

p density, kilogram mass/cubic meter 

E summation 

a standard deviation of Gaussian beam, centimeters 

T thickness ratio, dimensionless; also, vibrational 

relaxation time, seconds 

$ viscous dissipation function. Watts/cubic centimeter 

<p velocity potential, dimensionless 

viscous force vector. Newtons/cubic meter 
angular velocity, radians/second 

Superscripts 

( ) ' prime, dimensionless quantity 

( ) bar, normalized quantity; also average value of a 

quantity 

("') tilde, perturbation quantity 
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Subscripts 

a,b indicate values of a quantity in front of and 

behind a shock wave 

i x-direction or streamline index - first subscript - 

(i = 1 denotes freestream conditions) 

j y-direction or normal index - second subscript - 

(j =1 denotes centerline conditions) 

L linear value 

0 total or stagnation condition; also, peak value of 

an intensity distribution 

w wall or surface value 

00 freestream condition 

1,2, first, second, etc, perturbation quantities 

Miscellaneous 

V Del operator, 1/centimeters 
2 

V Laplace operator, 1/square centimeters 
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INTRODUCTION 



The optical quality of a propagating laser beam is de- 
graded when there are density gradients in the medium. A 
laser beam propagating through an absorbing medium creates 
local changes in the density which affect the index of refrac 
tion of the medium. Therefore, the density gradients will 
refract the light into regions of higher index of refraction 
(higher density) causing defocusing of the beam. This pro- 
cess, which is known as thermal blooming, has been studied 
by many authors [Refs. 1, 2, and 3] to determine its effect 
on a laser beam propagating through various media. 

When a laser beam is slued, a complicated coupled pro- 
cess of thermal expansion and fluid flow (relative winds) 

A 

takes place. Depending on the rate of sluing, the beam will 
experience subsonic, transonic, and supersonic winds at vari- 
ous radial distances from the laser source. The subsonic 
and supersonic thermal blooming problems have been 
solved [Refs. 4-8]. A computer program (DIDER) has been 
developed and experimentally verified [Ref. 9] for beam 
interaction with supersonic flows internal to a Gas Dynamic 
Laser. 

A solution for the transonic regime, and, in particular, 
the case when the frees tream Mach number is precisely unity, 
has been of considerable interest because the addition of 
heat can lead to extremely strong density gradients. As 
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with the theoretical investigation of transonic flows over 
airfoils and wedges, the linear small perturbation flow 
equations valid for subsonic flow (elliptic equations) and 
supersonic flow (hyperbolic equations) are no longer complete 
in the transonic regime. 

Many authors [Refs. 10 - 15] have studied various methods 
of solving the transonic flow equations without heat addition. 
The assumptions based on known transonic behavior about air- 
foils and wedges are used to simplify the governing equations 
and, further, to obtain approximate numerical solutions. 

Often, the assumptions greatly restrict the types of transonic 
flows that can be solved. For example, in the parabolic 
method [Ref. 16], the transonic flow equation is transformed 
into a parabolic differential equation. Since the solutions 
to parabolic differential equations are well established in 
the literature, this method may be used when the acceleration 
or deceleration is approximately constant in the region of 
interest (i.e. across the cord of many airfoils) ; however, 
this method is not valid in the region of the stagnation 
point or where the flow acceleration changes sign. 

With heat addition, the basic flow equations for transonic 
flow become more complicated. Zierep [Ref. 17] has solved 
the case when the Mach number is precisely unity; but, in 
solving the equation, he has reduced the problem to that of 
one dimension. Others [Refs. 18 and 19] have solved the 
nonsteady transonic equation with heat addition; however, 
the results are valid only for a short time after beam turn 
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on (i.e. no steady state solution) . Even for the short time 
solutions, extrapolation techniques between the property 
values at the critical points on the subsonic side and the 
supersonic side are employed to obtain the values across 
the transonic regime. Experimental investigations are being 
conducted [Ref. 20] to determine the influence of transonic 
sluing on thermal blooming. These experiments indicate that 
a steady two-dimensional solution is possible. Further 
experimentation will provide better guidelines for a 
theoretical analysis of the problem. 

In this thesis, transonic thermal blooming in steady 
two-dimensional invisid flow is formulated and solved for a 
Gaussian heat release distribution. In particular, it examines 
the problem when the freestream Mach number is precisely 
unity. This is accomplished by using a method developed by 
Broadbent [Refs. 21 - 25] in which the streamlines through- 
out the flow are adjusted in such a manner that the heat 
addition necessary to achieve these streamlines equals the 
required heat distribution and boundary conditions. The 
method is exact since specifying the streamlines reduces 
the nonlinear coupled momentum equations in natural (stream- 
line) coordinates to simple algebraic linear difference equa- 
tions that can be "marched" through the entire flow field. 
Boundary conditions are the freestream properties upstream 
of the heat release region and those on a boxonding stream- 
tube sufficiently removed from the heat release region that 
the flow is considered isentropic. The boundary conditions 
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on the bounding streamtube are calculated by methods similar 
to those employed in determining the velocity and pressure 
distributions over airfoils of known shape. Since stagna- 
tion points do not exist on the bounding streamtube flow, 
a modification to the method presented by Sprieter and 
Alksne [Ref. 10] for airfoils and wedges was developed. 

In Section II, the theory used in analyzing the transonic 
thermal blooming problem is presented and includes a des- 
cription of the assumptions made during the mathematical 
analysis, a summary of the pertinent equations derived in 
the Appendices, and a discussion of the method of solution. 

In Section III, the numerical results for the Mach 1.0 flow 
with a Gaussian heat release distribution are presented. 
Lastly, Section IV gives a brief discussion of the conclu- 
sions and of the other methods that show promise in solving 
the transonic thermal blooming problem. 
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II. THEORY 



A laser beam that is slued through the atmosphere at a 
constant rate of radians per second will experience regions 
of subsonic, transonic, and supersonic flow in planes per- 
pendicular to its axis at various radial positions along the 
beam. Fig. 1. The region to be investigated in this paper 
is the transonic region in which mixed subsonic and super- 
sonic flow exists, as in contrast to purely subsonic flow or 
supersonic flow, and, in particular, the case when the Mach 
number is precisely unity; that is, at a radial position 
given by Eq. (1) . 



/yRT, 

T — 7 ** = ■ 

M^ = 1.0 ft 



( 1 ) 



For this analysis, the laser beam is initially of a given 
axial symmetric Gaussian intensity sluing into an isotropic, 

-3 

quiescent medium of constant density p(kgm m ) . Viscosity, 
and thermal conduction are neglected. Justification for 
this assumption is given in Appendix A. The heating effect 
of the laser beam on the medium is approximated by a molecular 
relaxation process which is rapid enough that the absorbed 
energy is instantaneously transformed into heat. The energy 
distribution is of the form given in Eq. (2) . 



Q = la 



Iq a exp (- 



(x^+y 

2o^ 




(2) 
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FIGURE 1. Flow Regimes Along a Sluing Laser Beam 
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where Q (Watts cm ) is the time rate of heat absorbed per 

-2 

volume; Iq (Watts m ) is the beam maximum intensity; 
a (cm is the absorption coefficient; and the exponential 
of the Gaussian with standard deviation a (cm) , which is 
related to the spot size of the beam by w = xpl a, is 
centered at the origin of the x-y coordinate system trans- 
verse to the beam axis at r*. 

The theory employed for this investigation separates 
the steady transonic flow with heat addition into two en- 
tirely different flow regimes, Fig. 2. The internal, heat 
addition region can be treated exactly by the Broadbent 
method [Ref. 18] , in which the nonlinear equations of motion 
and the energy equation are solved numerically in natural 
coordinates. Then, a mesh of streamlines and normals is 
constructed throughout the flow field. Fig. 3., which reduces 
the momentum equations to simple algebraic equations solvable 
for the flow properties. Finally, the mesh of streamlines 
is adjusted until the properties calculated from the linear- 
ized momentum equations satisfy the energy equation and 
boundary conditions. Boundary conditions are obtained from 
the external, isentropic flov; region by approximating the 
transonic solution over a bounding streamtxabe to the internal 
region. 

The governing equations in the internal region have been 
derived for flow with heat addition, but without viscosity, 
heat conduction or body forces. A small perturbation analysis 
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FIGURE 2. Transonic Flow Configuration with Heat Addition 




FIGURE 3. Natural Coordinate System with Flow Mesh 
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for transonic flow (Appendix A) shows that viscosity, heat 
conduction, and body forces can be neglected to first order. 
The small perturbation, steady two-dimensional transonic 
equation with head addition is written as; 



(1-M 



^ 'f’x'x' 



+ <p 



y ' y ' ^co 






1 8q(x' ,y') 

CpT JiP 

(3) 



where the term on the right-hand side is added for flows 
with heat addition; q(x',y') (Joules kgm is the heating 
function; Cp (Joules kgm ^ °K is the heat capacity; 

T (°K) is the freestream static temperature; M is the 

CO CO 

freestream Mach number; and (j) is the nondimens ional velocity 
potential. Therefore, in vector notation, the steady state 
equations [Ref. 26] are; 



V (pu) =0 


(4) 


pU'VU + Vp = 0 


(5) 


pU*V(h + |u*U) = Q = la 


(6) 


p = pRT 


(7) 



— -1 -2 
where U (m sec ) is the velocity; p (newtons m ) is the 

static pressure; h (Joules kgm is the enthalpy; and 

R (Joules kgm ^ ®K) is the universal gas constant. Using 

a two-dimensional natural coordinate system [Ref. 26] , 

shown in Fig. 3, the governing equations were derived from 
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the above equations. The momentum equations are represented 
in the following convenient n\imerical form (Appendix B) . 



u: 



- u! 



i+l,j+l i,j+l 



An ' . . . 



YM„ 



i+l,j+l ^i,j+l 



-P- = 0 



(8a) 



P - 



YM^ 



-P,' 



■(u.',, .^,v! 



+ u,' , J =0 



i+l/j+1 ^i,j+l 2 i+l,j+l i+l,j+l i+l,j i+l,j 

(8b) 

where An.' . is the streamtube width at the i,jth station; 

Y is the ratio of specific heats; v! .is the inverse curva- 

ture defined as An,' ./R; and the prime (') denotes dimension- 

less quantities with respect to the freestream values. The 

quantities with a bar and parenthesized subscripts represent 

average values between the indicated stations (i.e. An). . . 

is the average streamtube width between station i,j and i+l,j), 

and the properties with the tilde {") denote perturbation 

values from freestream conditions (i.e. u. . = (u. ^ ” u, •)) . 

i/D i/J 

The natural coordinates s and n are tangent and normal to 
the streamlines, respectively. Eqs. (8a) and (8b) represent 
two nonlinear equations since u', p', An' and v' are dependent 
variables. However, specifying the streamlines of the flow, 
these equations become linear equations solvable for 
and 

at point d and uf . . , , at point b. Fig. 3. 

1,3 +1 1 , 3 J - 

The remaining flow properties and can 

be calculated from the continuity equation, Eq. (9) , and 
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equation of state, Eq. (10) , derived in Appendix B. How- 
ever, the new properties at the point e (i+l,j+l) may not 
satisfy the energy equation, Eq. (11) . If it is not satis- 
fied, an iteration process is performed by varying the 
streamline distribution and recalculating the properties 
until coincidence exists throughout the entire flow field. 
Hence, the solution of the problem is reduced to finding 
the streamline distribution that satisfies Eqs . (8) through 
(12) with the boundary conditions determined from the 
external region. 



Pi+1, j+l^i+1, j+A+l, j+1 



= 1 



Pi+l,j+l Pi+1, j+l^'^i+1, j + 1 



(9) 

(10) 






( 11 ) 



(laAn') 



(i,i+l),j As! . 



i, j+1 



where 









An 



l,j 



(Y“l) An^L^. 



(12) 



Insight into the behavior of the streamlines at the 
sonic speed was obtained from the solution of the linearized 
small perturbation equation with heat addition [Refs. 4- 
8] . A computer program (BLOOM) was devised (Appendix G) 
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which gave solutions for a flow very slightly below the sonic 
speed (M^ = 0.999), Figs. 4 and 6, and very slightly above 
the sonic speed (M^ = 1.001), Figs. 5 and 7. Computation 
of the nonlinear term in Eq. (3) from the linear results 
indicate that the linear equations are valid to = 0.9999 
subsonically and = 1.0001 supersonically. From the linear 
solutions, a trend of the streamline shape is obtained in 
the near sonic regime, even if not at sonic conditions. 

Nature usually exhibits a strong propensity to accomplish 
changes in a rather smooth fashion (even a shock wave is a 
smooth change in conditions if viewed microscopically 
enough); thus, a knowledge of flow conditions in the near 
sonic range gives a good clue to the sonic flow behavior. 

In Figs. 4 through 7, the Gaussian heat distribution is 
centered at the origin on the centerline, and the relative 
airflow is from left to right. 

From the appropriate density calculations and plots, 
discussed later in the thesis, it can be seen that the wake 
term becomes less significant compared with the source term, 
which is directly associated with the horizontal velocity 
perturbation, as sonic velocity is approached from either 
direction (i.e. M < 1 and M > 1) . Figs. 4 and 5 show the 
horizontal (longitudinal) velocity perturbation at Mach 
numbers of 0.999 and 1.001, respectively. Figs. 6 and 7 
show the vertical (transverse) velocity perturbation at these 
same Mach numbers. Note that, subsonically, the horizontal 
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FIGURE 4. Velocity Perturbation (u*) in Flow Direction for M^= 0.999 
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FIGURE 5. Velocity Perturbation (u') in Flow Direction for M^= 1.001 
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FIGURE 6. Lateral Velocity (v') Perturbation for M = 0.999 
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FIGURE 7. Lateral Velocity (v*) Perturbation for M = 1.001 
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perturbation. Fig. 4, decreases to a negative valley 
slightly upstream of the laser beam center and then changes 
rapidly to an equal positive peak slightly farther downstream 
of the beam center. The perturbation contours are of very 
small curvature near the beam center; the supersonic per- 
turbations, Fig. 5, show a negative valley near the beam 
center with a return to freestream conditions downstream of 
the center. In the supersonic case, the contours are nearly 
normal to the flow throughout, in keeping with the character 
of supersonic flow; whereas, in the subsonic case, the con- 
tours are closed, in keeping with subsonic flow behavior. 

The patterns of disturbance, however, appear to be consistent 
as the sonic speed is crossed. 

The vertical or transverse velocity perturbation. Figs. 

6 and 7, show similar behavior in the two speed regimes. 

As expected, the subsonic flow perturbation profiles are 
closed while those for supersonic flow are open to infinity. 
In both cases, the perturbations exhibit an odd symmetry 
about the main flow centerline. For the subsonic case, the 
perturbations are greater than zero to the left of the rela- 
tive flow and less than zero to the right. This is a direct 
indication of the spreading of the streamlines as the heat 
release region is approached by the air flow. Downstream 
of the laser beam, the streamlines again become parallel 
with the undisturbed flow direction. The supersonic case 
shows the same patterns except that the spreading effect on 
the streamlines is felt infinitely far in the transverse 
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direction, roughly along the characteristic lines from the 
disturbances . 

The effect of spreading of the streamlines as the beam 
is passed, which can be constructed from these plots of 
longitudinal and transverse velocity perturbations, is thus 
seen to be quite similar immediately below and immediately 
above the sonic speed. Fig. 2 is a sketch of the stream- 
lines. It is expected therefore that the streamline behavior 
at the sonic speed will not vary radically from this pattern. 

Based on the streamline distribution in Fig. 2, the 
following relationship between area change (streamtube width 
variation) and change in total temperature (heat input) was 
chosen, Eq. (13) . 



dA _ ^ 
A n 



~(y+l) (1+6) 




(13) 



where is the local total temperature of the flow and 6 
is a variable. One dimensional influence coefficients with 
area change and heat addition [Ref. 27], Eq. (14), were 
used to determine the required 6 at precisely Mach 1.0. 



dM 



2 



33T 



G(x) 

(1-M^) 



(14) 



where 



G(x) 




(1+|(Y-1)M^) (-2 



din (n) 
dx 



+ (1+yM^) 



dln(T^) 

dx 
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In order for Eg. (14) to remain finite and to be defined 
at Mach 1.0, G(x) must tend to zero in the limit as the Mach 
number approaches unity. If Eg. (13) is substituted into 
Eg. (14) and the limit taken, it is found that it must be 
zero at precisely Mach 1.0 (Appendix C) . Hence, the stream- 
tvibe area change is related to the heat input by Eg. (15) . 



dn! . 

n ' T 

1.3 



dT' 

l(,.l) ^ 

°l.j 



(15) 



where the maximum difference between n! . and n' ., and T' 

1.3 1.3 o. 



and T' 
o 



.-6 



1.3 



1.3 



is of the order of 10 and, hence, the freestream 



values can be used in the denominator. Using this approxi- 
mation for the streamtubes, an explicit formulation of the 

streamtube shape, at any position s! . along a streamtube, 

1.3 

was obtained (Appendix C) . 



" 1,1 ' 



(16a) 



1 ~ 

n! . = -yAn! , + E An! , + (j 

i,j 2 1,1 1,1+1 






i > 1. j > 2 



(16b) 



where for a Gaussian heat input at Mach 1.0 (Appendix C) , 




(Y+1) yia'I a 
o 



2a' 



s f 

1+erf ( — i-) 
J2a' 



i ^ 1, j 1 1 
(16c) 
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Once the streamtxabe shape is known, the slope (dy/dx) , 
Eq. (17) , and the radius of curvature (R) , Eq. (18) , can 
be calculated at every mesh point throughout the flow. 




1.3 




j > 1 



(17) 




*^i/3 ~^i-l,j 
As ! ! 



i > 1 



3 1 1 
(18) 



Therefore, the proportionality between the streamtxabe 
shape and heat release distribution gives not only an 
explicit formulation for the streamtubes and, hence, linear- 
ization of the nonlinear momentum equations but also simul- 
taneous agreement between those properties calculated from 
the momentum equations and the energy equation. The problem, 
in the case of precisely = 1.0 flow, becomes one of 
determining the correct boundary conditions . 

The boundary conditions that are used are those corres- 
ponding to freestream conditions (ui . = p' . = pi . = Ti . - 0) , 

1.3 1.3 1.3 1.3 

which are far enough upstream that the properties are not 
affected by the perturbation caused by the heat release 
region, and the pressure (p! .*) and/or velocity (u! .*) on 

the j*th streamtube which is considered far enough removed 
from the heat release region that the flow can be considered 
isentropic. Eqs. (16b) and (16c) can be used to determine 
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the streamtube shape and displacement where j = j*. However, 
a more convenient form is derived in Appendix C for a 

Gaussian heat release distribution at M =1.0 which can be 

00 

written as: 



n 



I 

i/ 




(Y + 1) Iq^C ' ^7T 

40 



s ! 

1+erf ( — 

y2a’ 



i ^ 1, j ^ 4 . 5 
(19) 



and the slope of the bounding streamtube: 



dn! .* 

. y ,;* = 

diJTjT 



(y+l) IqU yirfa ' 



,2 



4Q 



-exp - 



f 



2a 






exp - 



2a 



TT 



( 20 ) 



where t' is the maximum growth in the streamtube. 

Since the mass flow is constant in a streamtube, the 
bounding streamtube can be considered as a solid wall simi- 
lar to the surface of a wedge with a cusp extending forward 
to the freestream conditions upstream at infinity. Fig. 2. 
By hypothesis, this flow is isentropic and can be treated 
by methods similar to those used in analyzing transonic 
airfoil and wedge flows. In particular, an approximate 
solution is obtained through an iteration process on the 
integral equation derived by Green's f;anctions for small 
disturbance transonic flow similar to that performed by 
Sprieter and Alksne [Ref. 10 and 11] for thin airfoils and 
wedges but modified for the streamtube configuration which 
has no stagnation points. The absence of stagnation points 
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in the case of the flow around the bounding streamtube 
enhances the accuracy of the approximate numerical tech- 
niques because the stagnation point represents a singularity 
in the problem and places added constraints on the method and 
accuracy of the solution in the region of the stagnation 
point. In this method, the quadratic, nonlinear nature of 
the governing equations is retained which precludes shock- 
free supercritical flows in the transonic regime. 

The integral equation for transonic flow is derived by 
a Green's function analysis on the small perturbation 
transonic equation without heat addition, Eq. (3) without 
the right-hand side, resulting in Eq. (21) . 




(21) 



w 



where 




M ^(y+1)u' (x,0) 
00 * 



M 



00 



2 



( 22 ) 



2 



Ut (x,0) 



IT 



1 




(23) 



I 




(24) 




2it 



0 




( 24 ) 
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and the normalized variables are defined in Appendix D. 

Eq. (21) is a function only of x, the normalized coordinate 
along the wall surface, and t, the normalized amplitude of 
the wall slope, given by: 



-3 = T exp - (2 5a) 

dx 2 a^ 



where 

M^(Y+1)t M^^(y^- 1)I aJ2F a 

T = — = = ^ 2 (25b) 

3 PooV 

The iteration technique (Appendix D) is performed for 
various subsonic freestream Mach numbers starting at the 
initial supercritical flow and proceeding toward the Mach 
1.0 flow. Each solution obtained by the iteration technique 
gives a unique x , which determines the freestream Mach number 

and its corresponding normalized velocity distribution 

— — — --2/3 

u^(x,0) along the surface. If (u^-1) versus x is plotted 

for each position along the surface, it is found that as 

the freestream Mach number tends toward unity the slope 

becomes constant. This corresponds to the phenomenon of the 

Mach number freeze [Ref. 26] wherein the local Mach number 

is invariant with changes in the freestream Mach number 

when the latter is near unity or, more precisely. 



dM 

dM„ 

00 



= 0 



(26) 
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Vincenti and Wagoner [Ref. 28], Liepmann and Bryson [Refs. 
29 and 30] , and others have shown that the corresponding 
approximate relation yielded by the small-disturbance 
transonic theory is given by: 



d5 

^00 




(27) 



where 



5 



(1-M^) 



(u-1) 

7273 



(28) 



The freeze of Mach number extends over a finite range near 
Mach 1.0 where it is constant. Hence, the slope of the 
(u-1) versus t curve at each location x, is constant near 
the sonic condition, and the local Mach number at each 
position can be calculated by Eq. (28) . With the local 
Mach nxamber distribution known, the velocity perturbation 
can be calculated from Eq. (22) , and the other flow proper- 
ties can be determined from isentropic flow relations [Ref. 
27] . 

The Mach 1.0 flow is completely specified. Boundary 
conditions are known upstream of the heat release region 
and on a bounding streamtube. Calculation of the entire 
flow field is readily calculated from the simple algebraic 
equations obtained by Broadbent's method, where the stream- 
tube configuration is known precisely at Mach 1.0. 
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III. NUMERICAL RESULTS 



The computed flow properties are for the Mach 1.0 thermal 
blooming problem with the characteristics listed in Table I. 



Gaussian Beam 




Peak Intensity (I^) 


3.333xlo”^ Watts /m^ 


Standard Deviation (a) 


5 cm. 


Atmospheric Absorption (a) 


-7 -1 

8.0 X 10 cm 


Atmospheric Freestream 
Conditions 




Temperature (T^) 


288.0 «K 


Pressure (p ) 


1.013x10^ Newtons/m^ 


Density (p^) 


1.2246 kgm/m^ 


Velocity (u = a*) 

00 


340.0 m/sec 


Streamtube Width 

(An^^j=An^) 


1 cm. 


Streamline Cell Length 

(ASi_i) 


1 cm. 


Ratio of Specific Heats (y) 


1.4 


Specific Heat at Constant 
Pressure (Cp) 


1005.0 Joules/kgm^K 



Table I 

Freestream Flow Properties And 
Laser Beam Characteristics 
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The first step in the procedure is to determine the 
boundary conditions for the heat release region. Far up- 
stream of the heat release region, all the dimensionless 
perturbation quantities are zero. The other boundary con- 
ditions are calculated by Sprieter's method [Ref. 10] for 
the bounding streamtube, the shape of which depends upon 
the heat release distribution and is known, Eq. (16) . 

Using Sprieter's numerical iteration procedure on the 

transonic flow integral equation, Eq. (21) , the normalized 

velocity distribution along the bounding streamt\ibe was 

obtained for several frees tr earn Mach numbers between the 

initial supercritical flow and Mach 1.0. The results are 

graphed in Fig. 8. From these solutions, a (u-1) versus 
—2/3 

T curve. Fig. 9, was plotted for several positions along 
the streamtube. The results indicate that the slope (5) 
is constant for each position and, hence, the Mach number 
can be considered frozen. 

From the values of C(x), the local Mach nximber variation 
at a freestream Mach number of unity was calculated from 
Eq. (2'8) , since 5 is known and constant for the known 
physical streamtube shape, Eq. (30) , from Eq. (D6) . 

T = — - — = 1.386 X lO"^ (30) 

JYn" a 



where 



^ _ ( Y~l) Lairg' 



YP u 

' *^00 00 



= 1.737 X 10 



-5 



cm. 
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(a) Shock Wave at x = 9.5 
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FIGURE 8. High Subsonic Supercritical Flows on 
Bounding Streamtube 
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FIGURE 9. Mach Nxunber Freeze C = Constant for 
Various x Positions 
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is the maxiimam growth of the bounding streamtube . 

Once the local Mach number was determined, the velocity 
perturbation along the bounding streamtube was calculated 
from Eq. (22) and is plotted in Fig. 10. The remaining 
flow properties were calculated by the isentropic flow 
relations [Ref. 27] . 

It has been shown that the area change of the stream- 
tiibes in the heat addition region is directly proportional 
to the change in total temperature; that is, B = 0 in 
Eq. (13) . The energy equation is automatically satisfied 
by this selection for the area variation. Since the heat 
release distribution is known, all of the streamtubes are 
specified. From the streamtube shape, the radius of curva- 
ture of the streamlines were calculated at every point 
throughout the region, Eq. (18) . 

With the curvature known, the calculation of the flow 
properties in the heat region is reduced to the solution 
of two simple linear algebraic equations, Eqs. (8a) and 
(8b), the momentum equations in natural coordinates. The 
results of this "marching” technique are shown in Fig. 11 
where the density perturbations in the upper half of the 
symmetric laser beam are presented. This method has the 
distinct advantage that it can be solved on an HP 9830 
computer. 

Figs. 12 and 13 show the linear supersonic and subsonic 
density perturbation solutions calculated from the computer 
program (BLOOM) developed in Appendix G. The results 
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FIGURE 10. Velocity Perturbations on Bounding Streamtube at M 
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FIGURE 11. Density Perturbation (p‘) for 
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FIGURE 12. Density Perturbations (p') for = 0,999 
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FIGURE 13. Density Perturbation (p') for - 1.001 
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presented in Fig. 10 at Mach 1.0 are consistent with a 
continuous transition through the transonic regime. Although 
the solution calculated is for precisely Mach 1.0, it can 
be logically hypothesized that the flow behavior through 
the transonic regime is as shown graphically in Fig. 14. 

The linear supersonic and subsonic solutions , and 
Sprieter's high subsonic solutions are in excellent agree- 
ment with this hypothesis as well as that conjectured from 
the hodograph plane representation for transonic flow 
(Appendix E) . Therefore, there is a steady, two-dimensional 
solution for the transonic thermal blooming problem which 
is finite and does not result in "catastrophic" defocusing 
of the laser beam. 
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IV. CONCLUSIONS 



It has been demonstrated that there is a two-dimensional 
steady state solution to the transonic thermal blooming 
problem for a sluing laser beam with a frees tream Mach number 
of precisely unity. According to the linearized solutions 
for thermal blooming, the density perturbations become 
infinite as a Mach nximber of unity is approached. Due to 
nonlinear effects the trend to infinity is cutoff, and the 
finite value of perturbation quantities has been established. 

In obtaining the solution, it was found that at precisely 
Mach 1.0 that the streamtube area variation through the heat 
release region is directly proportional to the change in 
total temperature in that region. With the streamtube 
shapes known, an exact solution for the flow field can be 
obtained by Broadbent's method [Ref. 21] with the corres- 
ponding boundary values for Mach 1.0 isentropic flow far 
from the heat release region. There are no restrictions on 
the existence of shock waves nor on the symmetry or asymmetry 
of the heat release distribution, provided the bounding 
streamtube in the isentropic region and the boundary 
conditions on it can be calculated. 

Furthermore, in obtaining the solution at precisely 
Mach 1.0, it was found that shock waves exist throughout 
the transonic regime. Fig. 14. As the sonic condition is 
approached from subsonic speeds, a shock wave forms at the 
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FIGURE 14 . Laser Beam Centerline Mach 
Number Distribution in the 



downstream 1/e (one spot size) size, of the Gaussian heat 
release distribution, where the local Mach number first 
becomes sonic, and rapidly moves downstream to infinity. 

If the sonic speed is approached supersonically, a shock 
is formed at the beam center. To be consistent with the 
results in the subsonic portion of the transonic regime, 
it has been hypothesized that the shock moves rapidly up- 
stream to infinity. Finally, the analysis in Appendix F 
indicates that there is a minimum laser sluing rate that 
must be maintained to preclude significant phase distortion 
in the beam at Mach 1.0. 

The results of the Mach 1.0 thermal blooming problem 

for a Gaussian intensity distribution are in excellent 

agreement with the approximate results obtained by Ellinwood 

and Mirels [Ref. 19]. Fig. 15a shows these results where 

the maximum density perturbation was calculated to be of 
2/3 

the order of Q . Q is defined as: 

(Y-1) I ar 

Q = ^ , (31) 

Ypa 

where r is the radius (spot size) of a circular beam of 

uniform heating intensity (Iq^) . This approximation, when 

applied to the Gaussian beam of Section III, gives a maximum 

-4 

density perturbation of the order of 10 which is consistent 
with that obtained in this paper, and shown in Fig. 15b. 
Although the density perturbations are very small (i.e. 
measured in parts per million) , significant optical defocusing 
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(a) Results cf Sllinvrocd and Mirels [Ref. 16] 




(b) Results of Linear Sibsonic and Supersonic Solutions from DIEER 
[Refs. 4 and 5] and of Broadbents Method at Mach 1.0 

FIGURE 15. Schematic Representation of the Steady State 
Density Perturbations for a Sluing 
Two-Dimensional Laser Beam 
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of the laser beam could result (Appendix F) . With an 
increase in the beam area due to def causing, a large decrease 
in the beam intensity I (Watts/m ) will occur. 

Two other methods that show promise in solving the tran- 
sonic thermal blooming problem are the finite difference 
method developed by Murman and Cole [Ref. 12] for isentropic 
transonic flow which can be extended to include heat addi- 
tion, and the "finite element method" of solving partial 
differential equations [Ref. 31] which is a powerful method 
for solving a wide range of engineering problems . 

A Mach 1.0, steady state solution for thermal blooming 
in a laser beam with a given Gaussian heat release distribution 
has been presented in this thesis. Fig. 11. Subsequent inves- 
tigations should be performed to ascertain the Mach 1.0 
solutions for various Gaussian and non-Gaussian heat release 
distributions, beam intensities (I), atmospheric absorption 
coefficients (a) , and freestream flow properties. In this 
manner, a correlation might be developed that would be use- 
ful as an engineering approximation to the density perturba- 
tions at various transonic Mach numbers, vice a simple order 
of magnitude linear extrapolation as presented by Ellinwood 
and Mirels [Ref. 19] . Although solutions are known at Mach 
1.0 and in the linear subsonic and supersonic regions, from 
the computer program BLOOM, it would be premature to estab- 
lish a criterion that could be used in estimating transonic 
behavior in general before additional work is completed. 
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APPENDIX A 



DERIVATION OF STEADY TRANSONIC FLOW 
FOR A CONDUCTING GAS WITH HEAT ADDITION 

Eqs. (Al) through (A4) govern steady compressible flow 
with heat addition and constant properties [Ref. 32] . These 
equations include viscosity and heat conduction. The 
distribution of the heat addition is given by the function 
Q^f(x,y) and is assumed known. The gas flow is assumed 
iiniform at infinity (x = ±oo) . 



V* (pU) = 0 



(Al) 



pU‘VU + Vp = 1]; 



(A2) 



where 



= nV(v*u)/3 + nv^u 



pTU* Vs = ~ - V *q 



(A3) 



where 



^ (vu)^ - n(v^(U‘U) - (vxu)^ ~ 2u-v^u) 



and 




2 2 
Vp = pa Vs/Cp + a Vp 



(A4a) 



or 




(A4b) 
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The following dimensionless variables are defined; 



X' 


= x/^3^ 


y' = Y/I 2 




p' 


8 

Q. 

Q. 

II 


T' = T/T„ 


u' = u/u 


s * 


= s/Cv„ 


n' = n/n„ = 1 


k' = k/k 
' < 



where and I 2 characteristic lengths in the x and y 

directions, respectively, and the reference values are the 
freestream conditions. Substituting the dimensionless 
variables into Eqs. (Al) through (A4) and eliminating the 
pressure term in the momentum equation, Eq. (A2) , with the 
equation of ’State, Eq. (A4) , the equations simplify to the 
following; 



p* V •U' + U' 


. V ‘p ' = 0 




(A5) 


p'U' .V'U' = 


T' V p ' 


p'T'V's' V'(V'-U') V'^U' 


(A6) 


1 

CN 

1 


2 3 Re Re 

* 00 


p'T'U* ’V's* 




1 - "I (V -u') ^ + v'^(u' *u) 




Re 








- (VXU')^ - 2(U' 






YV'^T' 

RePr 


■ + Y(Y-l)M^^H„f (X,y) 


(A7) 



where the dimensionless ratios are defined as follows: 



Reynolds number: 


Re = p U t->/r\ 


Mach number; 


M = U //yRT 


Prandtl number: 


Pr = n /k 

n '00' 00 

*^oo 


Heat quantity: 




Scaling factor: 


L = ^,/^2 
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For flows in which conduction and viscosity can be 
neglected, two flows will be dynamically similar when y, 
(Mach number in the absence of heat addition) and 
are the same. 

Retaining the viscous, conduction, and heat addition 
terms, the governing equations, Eqs. (A5) through (A7), in 
two-dimensions become; 



p' (|^ + L 1^) + u' + Lv' = 0 



3x' 



3y 



(A8) 



x-direction momentum equation: 






2„. 



T' 3p' _ p'T* 3s' ^ _4_ 3^u 
~ 2 3x' 2 3x' 



3Re ^^,2 



3^v' 



3Re 3x'3y' Re , 2 



3^u' 



(A9) 



y-direction momentvim equation: 



It I 3v' , T I 3y^\ 

f aF- 3F-’ 



energy equation; 



p.T'(u' + Lv h;) 



LT' 3p' 


Lp'T' 

ym ^ 


3y' ^ 


4L^ 3^v' 

3y.2 


L 3^u' 


+ 1 


3^v' 


(AlO) 


3Re 3x'3y' 


' Re 


3x'^ 


y(y-i)m^^ 

Re 


1 


,3u*. 2 
3x ' ^ 


.. 3u'3V 
^ 3x'3y' 



3^V 



3y 



,2 






3u ' ^ 2 
3y ' 



+ Ri^ ^^772 ^ ^ 



3x' 



Y(Y-1)M^ f(x,y) 



3y' 
(All) 
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Since the transonic flow region is of interest, a 
transonic expansion procedure similar to that used for 
incompressible flow around an airfoil of thickness ratio t 
[Ref. 33] will be used. In the following analysis, the 
heat addition will cause the perturbations from the free- 
stream conditions. Introducing the following perturbation 
quantities : 



u' 


= 1 + 


2/3 ^ 4/3 

e u^ + e ' U 2 ; 


v' = 


ev^ 


p’ 


= 1 + 


2/3 ^ 4/3 

e p| + e ' p^ ; 


T' = 


1 + e^/^ + e'^/^ 


s ' 


= 1 + 


§2 ; 1/Re ~ 0(e) 


f 


H„~0(e^/^) ; 


"^1 


= 1 


; I 2 = ; L 


= 




the 


first and second order equations 


at M„ 

00 


= 1.0 become; 




First order equations : 








9u| 

aF" 


ap,’ 

^ 93T= ° 




(A12) 




3vi 

3x' 


3p^ 

+ = 0 
ay' 




(A13) 




9T| 


9uJ 

+ (Y-l) ^ = 0 




(A14) 



Integrating Eq. (A12) gives u^ = which, when substituted 

into Eq. (A13) , illustrates that for first order theory the 
flow is irrotational . 
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0 



(A15) 



3v| 3u| 

aP" “ dp" 



Finally, Eq. (Al4) gives the temperature distribution as: 

T| = (Y-1) = -(Y-1) u{ (A16) 

Second order equations: 



3uA 3uJ 


3vi ; 


- I— 1 

Q. 


3y' : 


3u- ^ 3u- 


3u' 


’^l 3x^ 


Pi ^ 



3pA 3p^ 

^-4 + u' = 0 

3x' 1 3x' 



3P^ 



dp{ 



+ T' 



M. 



2 '3x 



2^ 3s 2 4e ' 3 uj^ 

rr =3 



ym 



2 3x' 



3x 



,2 



3v' 3v' 

Pi ^ ^ 



M 



. 3 Pi 

^^i 3y^ 



3p'^ 



ym. 



j_ Hi . 1/3 



2 3y' 



3x 



.2 



3|. 



3x 



r = Y£ 



1/3 



2^. 






3x 



4 + Y(Y-l) m 4 f (x,y) 

I ^ 



(A17) 



(A18) 



(A19) 



(A20) 



Combining Eqs . (A17) through (Al9) and using Eqs . (A12) and 

(A16 ) f the transonic perturbation equation is obtained. 

)^T' 

3x' 



4e^/^ 1/3 ^^^i 

-(Y+1) u| ^ 5 + e ^^^2 ■*■ (Y“l)f(x,y) 



(A21) 
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The first term on the right-hand side of Eq. (A21) 
represents the effect of viscosity, the second term the 



conduction, and the last the heat addition. As can be seen. 



are small and can be neglected under most circumstances. 

-2/3 

If, however, the Reynolds number is of the order of e 
or smaller, both terms become important and must be retained 
in the perturbation equation. Since to first order the flow 
is irrotational , the velocity potentials u| = 9<))/9x' and 
v| = 9(})/9y' can be introduced into Eq. (A21) giving: 



A similar analysis can be carried out for a freestream 
Mach number close to unity. In this case, the freestream 



. 1/3 

both the VISCOUS and conduction terms contain a e term 

when the Reynolds n\imber is of the order of e and, hence 



- ( Y+1) I I + 4>y I y I 



(3Y+l)e^'^^ ^ 

3 ^x'x'x' 



+ (Y-1) f(x',e 



(A22) 



2/3 

Mach nximber can be approximated by M^ = 1 + Ke ' and the 
transonic expansion gives: 







(Y-l) f (x' 



(A23) 
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APPENDIX B 



BRIEF DEVELOPMENT OF BROADBENT'S METHOD OF 
SOLVING NONLINEAR EQUATIONS' OF MOTION WITH 

HEAT ADDITION 

Broadbent's method [Ref. 24] is an inverse solution in 
which the streamlines, normals, and boundary conditions 
along at least one streamline and normal are specified and 
the resulting flow field calculated. 

The governing equations for this derivation are as 
follows : 



V (pU) = 0 


(Bl) 


pU'VlT + Vp = 0 


(B2) 


pu.v(h + jU-U) = Q = la 


(B3) 


p = pRT 


(B4) 



In natural coordinates [Ref. 26], these equations 
become : 



(pu4n) = 0 


(B5) 


If * H = ° 


(B6) 


2 






(B7) 


pu (h + = la 


(B8) 
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A mesh is constructed over the flow field using stream- 
lines and normals which are chosen in advance as in Fig. 3. 
The cells specify the streamtube boundaries and the mesh 
lines the streamline slope and curvature throughout the 
flow field. With the mesh specified, Eqs . (B5) through (B7) 

are nondimensionalized using uniform freestream flow 

properties (u, •/ Pi T. ■), a constant uniform 

■^fj -*-/J -*-fj J-tj 

streamtube width upstream of the perturbations (An,' .) and 

-*- / 3 

constant Cp. 



n! .u! .An! . = 1 
Px,3 1,3^ x,3 



(B9) 



u 



- u! 



i+l,j+l i,j+l 



^^(i,i+l) .j , 5 , _ 5 , 1 

2 ^i+l,j+l ^i,j+l 






(BlOa) 



^i+l,j+l ^i+l,j 



- i YM„2(u! 



V ! 



i+l,j+l i+l,j+l 



+ u ! , , ■ V ! , , . ) 

x+1,3 x+1,3' 



(BlOb) 






(laAn ^ i+i) i 



p! . = p! . RT! . 

^x,3 ^x,3 x,3 



(B12) 



where is a freestream heat quantity defined in Eq. (12); 

Q! . is a dimensionless heating term; and v! . is a 
1 , 3+1 ^ x ,3 

dimensionless inverse radius of curvature. 
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The slope (dy/dx) and the radius of curvature (R) of 

the streamlines at a point b (i+l,j) are shown geometrically 

in Figs. 3. Due to the symmetry of the problem, the slope 

and the curvature of the centerline streamline are zero 

(i.e. for all stations i ^ 1 and j = 1) . Assuming that 

Ay. . = An. . and Ax. . = As . . (i.e. to first order the 
1/3 1/3 1/3 1/3 

mesh is rectangular) , and that the streamline shape is 
known explicitly, the slope in finite difference form is 
expressed as: 



dx 



1/3 



■ *^i/j 



n.' T .-l/2An!^T . 

1+1,] ^ 1+1/3 



As! T 
1/3 



n! . + 1/2 An! . 
1/3 1/3 



(Bl3a) 



or : 



'i.: 



^^i+1, 1 j^Jj^^i+1, 1 ■ 2^^i+l, j ■ 2^i, 1 



3-1 .. 1 

E An! . +^n! . 



i/j 2 i,j 



As! . 
1/3 



(Bl3b) 



where 






+ 1 ) 



(Bl4a) 



n: 



= -iAnl 



j-1 



i/j 2 i,l 



■ E An! , — oAn! , + 
1=1 



— TrAn! . 
2 1,1 






i > 1, (Bl4b) 
j >2 



The inverse of the radius of curvature (v) is defined in 
Eq. (B15) as: 
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I 

2 X 0 



I 



* ' i!»dw 



V 



(B15) 



d^y/dx^ 



^ (l+(dy/dx)V^ 



By nondimens ionalizing and substituting the slope obtained 
in Eqs . (B13a) or (Bl3b) into Eq. (B15) , the following 
finite difference equation is obtained; 



An 



^^i^ j ~ ^i-1^ 



V : 



As ! T 



1/3 



7l + (6. 

^ / J 



(B16) 



Because the magnitude of the slope 6 . .is never greater 
“ 6 

than 10 , it can be neglected in the denominator simplifying 

the expression for the curvature. Substituting Eq. (B13b) 
into Eq. (B16) and retaining only terms of first order, the 
curvature becomes 



v: . = 
1/3 



j-1 ^ I ^ . 3-1 .. 

,+ Z An.' 1 , --iAnl^T .-An! ,-2 Z An? , 
2 xtl,l itl,l 2 itl,j i/1 ^ i/1 






(B17) 



1 . 3-1 .. 1 . 

+ An! . +^An! . , + Z An! . , --^An! , . 

1,3 2 1 - 1,1 1 - 1,1 2 1 - 1,3 



/(As! .) 

-*-/ J 



If the mesh is not specified, Eqs. (B9) through (B12) 
represent five equations in five unknowns and are nonlinear 
since An', p' and u' are dependent functions. When the 
streamt^abe shapes are specified. An' and the curvature v' 
are known quantities and the equations become linear in p ' 
and u' . Therefore, Eqs. (BlOa) and (BlOb) represent two 
algebraic finite difference equations solvable for p' and u' 
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throughout the flow field, provided suitable boundary condi- 
tions are known. For example, if the values of p' and u' 
are known at points a (i, j) , b (i+1, j) and d (i, j+1) in Fig. 

3, the pressure and velocity perturbations can be calculated 
algebraically at point e (i+1, j+1). This procedure is con- 
tinued ("marched”) throughout the mesh until all flow pro- 
perties have been determined. 

The remaining flow properties, density, and temperature 
are then calculated from Eqs. (B9) and (B12) , the continuity 
equation and equation of state. Once this has been accom- 
plished, the energy equation, Eq. (Bll) , is checked to see 
if the properties calculated from the specified streamline 
shapes agree with the desired heat distribution. If it does 
not coincide, an iteration procedure is performed in which the 
streamline shapes are varied until agreement is obtained. 

Boundary conditions are needed along at least one stream- 
line (usually a wall) and one normal (usually freestream 
conditions) . In the transonic thermal blooming problem, the 
boundary values are chosen as the unperturbed freestream 
properties upstream of the heat release region and the 
velocity or pressure on a bounding streamtube which is far 
enough away from the heat release region that it can be 
considered isentropic. Fig. 2. With the pressure and/or 
velocity known on two of the boundaries, the solution for 
all the flow properties can be calculated as mentioned pre- 
viously. If a coundary condition is known at other points 
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in the flow field, the resulting properties from the 
"marching” procedure must agree with it or the specified 
streamtube configuration must be altered until coincidence 
between the boundary conditions and heat distribution is 
obtained throughout the flow field. 
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APPENDIX C • 



DERIVATION FOR THE RELATION BETWEEN AREA 
VARIATION AND HEAT ADDITION AT PRECISELY 
MACH 1.0 FLOW 



1. AREA CHANGE WITH HEAT ADDITION AT THE SONIC POINT 

From the influence coefficients for constant specific 
heat and molecular weight [Ref. 27] , the change in Mach 
number with heat addition and area variation is given by: 



dM^ = G(x) 

dx . 

1-M 



(Cla) 



where 



o 1 -y y d(lnT ) H M r> a ^ 

G(x) = M^(1 +j(y- 1)M^) [ (1-fyM^) ■ - 2^^=^], (Clb) 

d(lnT^)/dx and d(lnA)/dx are assvimed to be functions of x, 
the streamwise coordinate, and dx is always positive. Based 
on the relation between streamline shape and total temperature 
change obtained from the linear subsonic flow and supersonic 
flow. Figs. 4 through 7, the following relationship between 
area variation and total temperature change (heat addition) 
is chosen: 

= |(Y+1) (l+6(x)) ^ . (C2) 
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3 is a variable providing an additional degree of freedom 
in determining the relationship at Mach 1.0. Substituting 
Eq. (C2) into Eq. (Clb ) , G(x) simplifies to: 

o 1 00 ^ ( In T ) 

G(x) = M^(1 +j(y- 1)M^) ((M^-1)Y - (Y+l)e) ^ (C3) 



Since the denominator of Eq. (Cla) tends to zero in 
the limit as the Mach number approaches unity, the numerator 
must go to zero to permit a continuous passage from subsonic 
to supersonic speeds (or vice versa) . The condition, when 
G(x) passes through zero simultaneously with the Mach number 
becoming unity, is designated the critical point G(x) and, 
from Eq. (C3) , is equal to: 



lim G(x) E G*(x) = -J(y+1)^333^ 0 (C4) 

M -^-1 



Hence, for heat addition -d(ln T^)/dx > 0, 3 must equal 
zero or, equivalently, the area change is directly propor- 
tional to the heat addition (i.e. change in total tempera- 
ture) for steady flow from Mach 1.0, Eq. (C5) is: 



= i(V-l) 



ddn T^) 
dx 



(C5) 



2 

At the critical point, dM /dx is of the indeterminate 

2 

form 0/0. Therefore, the limiting value of dM /dx at the 
critical point is found by applying L 'Hospital's Rule to the 
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right-hand side of Eg. (Cla) . 



= (^)* = j a g /dx) * 

dx dx -(dM^/dx)* 



(C6a) 



or, equivalently. 



(M_) * = ^ 

^ dx ^ 




(C6b) 



2 

Eg. (C6b) suggests that (dM /dx) * can be positive or negative 
and the flow can become either subsonic or supersonic from 
the sonic point. The requirements on 3 follow from Eg. (C6b) 
where G(x) is given by Eg. (C3) . 



( 



dx^ 



lim 

M-»-l 

3=0 



dx 



|(Y+1) 



-(Y+l) (§|)* 






d(ln T^) 
dx 



(C7) 



Thus, at the critical point, the flow may be continuous 
through the sonic speed provided (dG/dx) * is negative and 
is not zero. For heat addition, Eg. (C7) can be solved for 
the critical value of (d3/dx)*. Eg. (C8) ; 



(M)* 



Y 



(dMI) 



Y + 1 dx 



(C8) 



and is graphed in Fig. 1C. 
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FIGURE 1C. B Variation versus Mach Number Variation 
for Sonic Flow 

From the sonic flow condition, S must be positive for 
supersonic flow and can be either positive or negative 
for subsonic flow. The change in 6 with Mach number squared 
must be greater than y/(Y+D for either flow condition, 
supersonic or subsonic. 

In conclusion, for sonic flow the area variation is 
directly proportional to the change in total temperature, 

Eq. (C5) , where the change of 3 from the sonic point must 
satisfy Eq. (C8) . 
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2. STREAMLINE CONFIGURATION FOR MACH 1.0 FLOW WITH A 
GAUSSIAN HEAT INPUT DISTRIBUTION 

In the first part of this Appendix the streamtube area 
change was shown to be directly proportional to the change 
in total temperature at precisely M^ = 1.0; see Eq. (C5) . 

To derive an expression for the change in total temperature, 
a slightly different form of the energy equation is used, 

Eq. (C9) , 






(C9) 



in terms of the total temperature T^. Assuming a constant 
Cp and a Gaussian heat intensity distribution in the form 
given in Eq. (2), Eq. (C9) becomes: 
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(CIO) 



where n| ^ represents the total displacement of the jth 
streamtube from the centerline (j=l) , and defined by 

Eq. (12), is evaluated at M^ = 1.0. 

Equation (CIO) is now substituted into Eq. (C9) to give 
an expression for a change in total streamtxobe growth over 
a streamline distance As! at M = 1.0. 
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or more conveniently for computational purposes; 
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(C12b) 



for streamtiibes off of the centerline (j > 1) . Equations 

(C12a) and (C12b) can be integrated to give an expression 

for the total streamtube perturbation from its initial 

width n ' . as : 
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(C13b) 



An alternate form of Eqs. (C13a) and '(C13b) which can be 
used in Eqs. (B13a) and (B13b) to simplify the expressions 
for the slope and the radius of curvature of the streamlines 
is as follows: 
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i^l, j^l (C14) 



where it has been assiimed throughout that the initial 

streamtube widths are uniform; that is, An,' • = 1 or 

J- / j 

An^ • = 0 for all j . 

An expression for a streamtube, which is far enough from 
the heat release region that it can be considered isentropic, 
is now derived. Carrying out the integration over y' in 
Eq. (C12a) and assuming that n! . > 40' an expression for 

1 f D 

the slope of the bounding streamtube is obtained; 



dn! . * 

... = 

ds . 



1/ j* 



( Y+1) I^a/Tif o' 

40:: 



exp - 






2o 



.2 



i>l, j>l. 



(CIS) 



72 



i im$ts 



Integrating Eq. (C15) with respect to ds| gives an explicit 
equation for the boxanding streamtube perturbation from its 
initial position n£ j* as; 



n! . * ' 
1/3* 


1 . 




(Y+1)I aa'^ s! 

= 45 (1 + erf p^) . (C16) 
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APPENDIX D 



INTEGRAL EQUATION DERIVATION FOR TREATING 
SMALL DISTURBANCE TRANSONIC FLOW WITH SHOCKS 
FOR AN INFINITE WALL WITH GAUSSIAN SLOPE 



1. GENERAL DERIVATION OF INTEGRAL EQUATION 

In this appendix, an integral equation is derived for 
the small disturbance transonic flow with shocks on an 
infinite wall with Gaussian slope. The small perturbation, 
steady, two-dimensional transonic equation without heat 
addition can be written as: 



(Dl) 

where the local Mach number is approximated by the freestream 



Equation (Dl) is valid only in regions where the necessary 
derivatives exist and are continuous. Therefore, the integral 
equation derivation includes an equation for the transition 
through a shock; the equation is the classical relation for 
the shock polar [Ref. 34], Eq. (D2) . 
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Mach ntomber plus a perturbation quantity (M (y+1)<|)„i). 
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where subscripts a and b refer to conditions ahead of and 
behind the shock. If a small perturbation analysis is 
carried out on Eq. (D2) similar to that performed in 
obtaining Eq. (Dl) , the following relation is found for 
the perturbation velocity components across the shock wave: 

(1-M^) (u^ - u^)2 + (v^ - VjI,) ^ (Y+1) (u^ + u^) (u^ - u^) ^/2 

(D3) 

Equation (D3) corresponds to the shock polar curve for 
shock waves of small strength that are inclined at any angle 
between that of normal shock waves and that of the Mach lines. 

The applicable boundary conditions for the integral 
equation analysis are that the perturbation velocity must 
vanish far upstream of the disturbance, Eq. (D4) : 




X=-oo 



3y* 




0 



(D4) 



and that the flow must be tangent to the wall surface, 
Eq. (D5) ; 



34 > 



3Z* 
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'=0 



®i, j* 



T exp - 7 

20 '^ 



(D5) 



where 9z'/9x' is the slope of the wall surface in the x' (s') 
direction which can be approximated by Eq. (CIS) ; and t is an 
equivalent thickness ratio for the surface defined as the 
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maximum wall growth t (cm) in the y-direction divided by a 

characteristic length /2rr o for M = 1.0. 
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Also, it is necessary to prescribe that the direct influence 
of a disturbance in the supersonic region proceeds only in 
the downstream direction. 

Sprieter and Alksne [Ref. 10] derive in detail the 
integral equation by a Green's function analysis for transonic 
flow around bodies having subsonic freestream velocities, and 
a local shock discontinuity. The nonlinear term in the 
differential equation was maintained. This derivation can by 
extended to the streamtube configuration of the boundary 
condition for the internal region. For the sake of brevity, 
only the important equations and pertinent assumptions will 
be given in the following analysis. 

Since the object of the analysis is to determine the 
pressure perturbation and/or velocity perturbation on the 
bounding streamtube, an equivalent equation for the velocity 
perturbation can be readily obtained by differentiating 
Eq. (Dl) with respect to x'; it is: 




2 



(D7) 
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For the Green's function analysis, it is advantageous to 
normalize Eq. (D7) with the following linear transformation; 



X = x' ; y = 3y' 

_ _ (Y+l)u’ _ M^^(y+1)u' 

u = 5 ; V = = (D8) 

6 3-^ 

where 3 = /l - 



In this way, Eq. (D7) reduces to the following: 



8^u ^ afu 

8x^ 8y^ 



V^u 




(D9) 



where u can further be defined by the approximation to the 

2 2 

local Mach number, Eq. (Dl) , namely M = (1 + (y+Du), as; 

M ^(v+l)u' - M ^ 

mmm 00 ' 00 

u = = =• (DIO) 

1-M^ 1-M^ 

00 oo 

From Eq. (DIO) , it is clear that for a subsonic freestream 
Mach number: u < 1 when the local flow is subsonic, u = 1 
when it is sonic, and u > 1 when it is supersonic. 

The Green's function analysis is a Neuman problem for 
finding a function u in a bounded region R, bounded by a 
finite number of closed smooth contours C, such that Eq. (D9) 
is satisfied in R, u and its first partial derivative are 
continuous in R (shocks are excluded by judicious selection 
of the boionding contour) and the normal derivative 3u/8y 
is prescribed on the boundary of R, Fig. ID. 



77 



iMoor 




r ^ , 




FIGURE ID. Integration Region for Green's Function Analysis 
on Flows Past Bounding Streamtubes at 
Transonic Speeds 



Green's theorem states that: 



/ (u7G‘n-GVu‘n) dC = // (uV^G - GV^u) dR (Dll) 

C R 



where the directional derivative on the left side are 
taken along the normal, drawn inward, to the curve C. 
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The Green s function G is constructed so that it satisfies 
the Laplace equation v^G = 6 (C-x^n-y+Yo) with ac/ 9 n = 0 on 

1 r, 

G(c,n/X,y-y^) = ^ ^ 

where 

/(x-c)2 + (y-y -n)^ 

R = 'f • - ° 

J{-x.-X)^ + (y+yQ+rF) ^ 

and ^ and n are the variables of integration while x and y 
are the coordinates of a point P. With the Green's function 
of Eq. (D12), Eq. (Dll) becomes: 



CO 00 2 —2 

u(x,y) - J J G (^) dn’d'^ 

-coO 3 5 



J (G - u 



dj - j (G ^ - u -^) 



3n 



3n 



dn 



T=o 



9n 



3n 



X=X 






(G is . s 3G 



8 n 



8 n 



dn + / (G - u -^) dc (D13) 



x=x. 



3n 



3n 



where the respective integration regions and contours are 
designated in Fig. ID, Integrating the remaining surface 
integral by parts twice with respect to " 5 , applying the 
Green's function boundary condition on (tT=0 or y=y^) , 
allowing the radius r of the curve C 2 go to infinity, and 
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decomposing the line integral over the shock wave into 
components parallel to the axes of the coordinate system 
yields Eg. (D14) . 



u(x,y) - j (x,y) = 



u rfu^^dCdn 



-M 



G 1- _ ^-_l-2j ^ 

. 9C ^ ^ 3l 



G— (u-iu^) - ^ 

. 3T ^ ^ al 



G 9u _ - 3G 
_ BtT StTJ a 



G - SH 
_ 3n 3n _ 



cos ( 



cos (n,5) 



n>n) 1 

n .T\ J 



dn 



(D14) 

where represents the linear wall velocity u when the 

wall shape is specified. 



uL = / G IH . J M 

^ -» 3n 9n 



dX 



n=o 



(D15) 



Sprieter and Alksne [Ref. 10] in their derivation delineate 
the advantages for integrating the surface integral over R 
by parts as was done in arriving at Eg. (D14) . The advantages 
are important enough that they are reiterated here. First, 
the double integral of Eg. (D13) shows a very strong influence 
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of the velocities in the region immediately surrounding the 

point P since they are multiplied by ln(R) . This influence 

is"- largely nullified in the double integral of Eq. (D14) 

because part of the region has a negative influence and part 

has a positive influence. The predominant influence in the 
. . 1 —2 

latter case is furnished by the term j u standing outside 

the integral. Next, the contribution of distant regions 

is diminished in importance in the double integral of 

Eq. (D14) since their influence varies inversely with the 

square of the distance, rather than with the logarithm in 

Eq. (Dll) . A further advantage is that the value of the 

double integral of Eq. (D14) is continuous through a shock 

wave rather than discontinuous as is the case with Eq. (D13) . 

Lastly, a point of great importance in the approximate 

solution arises from the fact that the integration by parts 

1 —2 

provides extra terms (-those containing y u ) in the 
integrals along the shock surface S which combine with those 
already present in such a manner that the contribution of 
these integrals becomes very small when the shock wave 
approaches a normal wave, as is usually the case at high 
subsonic speeds. 

Therefore, u must satisfy Eq. (D14) where Uj^ is given 
by Eq. (D15) and the shock relationship Eq. (D16) , which is 
the normalized form of Eq. (D3) . 

■ “b> ^ ■ ^b> ^ =1 "b> <“a • “b> ^ 
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In the case of a specified symmetric wall shape, Eq, (D14) 
can be considerably simplified as it permits the introduction 
of the boundary conditions which are specified at the outset 
of the problem, Eqs. (Dl7a) and (D17b) , into the integral 
over the wall surface. 



u(x,0) = Uj^^ = KNOWN 



and 



9u _ ^ 
3n 3T 






(Dl7a) 



(Dl7b) 



where Z represents the reduced ordinate of the surface, and 
T the reduced thickness ratio given by: 

(Y+1) T 

T = = (D18) 

Substituting the relationships of Eqs. (Dl7a) and (Dl7b) into 
Eq. (Dl5) completely determines at the outset of the analysis 
the linear wall velocity distribution: 



Ul„(x,0) 
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2. SIMPLIFICATION AND APPROXIMATE SOLUTION OF THE 

INTEGRAL EQUATION FOR TRANSONIC FLOW 

To obtain a solution to Eq. (D14) , some further 
ass\imptions must be made. They are; (a) all shock waves 
lie in a plane transverse to the beam, and (b) the shock 
waves are normal (i.e. normal to the local flow direction). 

These two assumptions reduce Eq. (D16) to an expression, 

Eq. (D20) , which permits an advantageous introduction into 
Eq. (D14) to eliminate the integral over the shock S, Eq. (D21) . 










(u-| u2)b 



(D20) 



0000 _2 2 r 

u(x,y) = Ub(x,y) + j u^ (x,y) + ^ ^ 

2 (J321) 



Equation (D21) will be used as the basis for obtaining 
the velocity distribution on the bounding streamtube. 
Approximate solutions of this equation are obtained numerically 
by starting with an assumed velocity distribution and 
iterating until convergence is obtained. If mixed flow 
exists, the initial velocity distribution must include a 
proper discontinuity complying with the shock relations of 
Eq. (D20) . 
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Due to the double integral in Eq. (D21) , an iteration 
process would be cumbersome unless it could be reduced to 
a single integral by introducing suitable approximations. 
Oswatitsch [Refs. 35 and 36] has shown that approximate 
knowledge of the velocity distribution in the vicinity of 
the surface can be expressed in terras of the local coor- 
dinate y, the ordinates of the wall shape F(x), and the 
desired but unknown velocity distribution u^(x, 0) on the 
surface. This permits the reduction of the double integral 
to a single integral. 

Oswatitsch has considered such an approximate relation 

in which the velocity u(x/y) starts from the value u^(x,0) 

at the surface with an initial rate of change given by the 

irrotationality condition, Eq. (D22) , and vanishing at large 

—2 

distances as 1/y . 
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The resulting relation takes the following form: 



u(x,y) = 



u^(x,0) 



(1 + (y-YQ)/b) 



(D22) 



(D23) 



where b is a function of x satisfying the irrotationality 

condition at y = y . Solving for b using the boundary 

o 

conditions of Eqs. (D17a) and (D17b) yields: 
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(D24) 



2u (x,0) 
w ' 

— 112 “ • 
d^Z/dx 

Substituting Eq. (D23) into the double integral of Eq. (D21 ) , 
integrating with respect to tT and setting y = results in 
a simplified integral equation for calculating u^(X/0). 

For computational purposes, the integral equation is expressed 
as : 



u^(x,0) 






(D25) 



where 



I = / dC (D26) 

— 00 

and the function E is : 

E(-i^) = E(X) =-F[flX| (5 - lOX^ - 

b (1+X^)^ ^ 

(1 - lOX^ + 5x"^)lnlXl -3^(1+X^) (25-71X^-X^-X^) ] 



Although I is a function of u^ and is unknown, it is informative 
to rewrite Eq. (D25) by solving for u^ in terms of I and la^, 
thus ; 

u^ = 1 ±Jl - (2u^^-l) = 1 ±Jl - L , (D27) 
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where 



L = 2Ut - 1 
Lw 

a known function since u^^^ is known once the wall shape is 
specified. It should be noted from Eq. (D27) that the 
discriminant must always be positive in order to obtain real 
values for u , thus I > L. Furthermore/ the choice of the 

>Af 

plus or minus sign determines whether the local velocities 
are subsonic or supersonic. A change in sign at a point 
where the radical is zero corresponds to a smooth transition 
through the sonic speed. A change in sign at a point where 
the radical is not zero corresponds to a discontinuous jump 
in velocity, namely, a normal shock from supersonic to sub- 
sonic flow when progressing in the flow direction. A dis- 
continuity in the reverse direction is inadmissable since 
it corresponds to an expansion shock, an impossible phenomenon 
which violates the second law of thermodynamics. 

3. COMPUTATIONAL TECHNIQUE 

In this section, a brief description of the computational 
technique used to solve Eq. (D25) will be presented. Sprieter 
and Alksne [Ref. 10] describes in great detail how a judicious 
selection of the initial velocity distribution, depending 
upon the surface shape and flow configuration (i.e. entirely 
subsonic flow or mixed flow with possible shocks) , is essen- 
tial in providing a rapid convergence to the solution. It 
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will be assumed for the following discussion, therefore, 
that a judiciously selected velocity distribution is known. 
The method of solution is as follows: 

a. From the known surface shape in normalized coor- 
dinates Z, calculate its first two derivatives with respect 

to the normalized flow coordinate x; that is, d^/dx and 
2 — —2 

d Z/dx . Since the normalization depends upon the freestream 
Mach number, which is unknown, the value of T, Eq. (D18) , 
will be a parameter in the iteration technique. 

b. Calculate the linear wall velocity u^^(x, 0) from 
Eq. (D19) . As '? does not depend on Xr it can be taken 
outside of the integral of Eq. (D19) giving an expression 
directly proportional to the parameter t. 

" 7 (D28) 

where (x, 0) is a function of only the position x and is 
fixed once the surface shape is specified. The quantity 
L follows immediately from Eq. (D27) . 

c. Assume an initial value of t consistent with the 
velocity distribution and calculate b(x) from Eq. (D24) . 

d. Since b(x) is now known and an initial wall velocity 

u (x,0) has been selected, Eq. (D26) can be integrated 
w 

numerically to get I . 

e. The I function calculated for the initial value of 
T is graphed, t is then varied in the L function until 
the two curves are tangent at, at least, one point (i.e. 
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I-L = 0) . As was stated earlier, I ^ L at all points x 
along the surface for the velocity to be real. 

f. Using the L function that satisfies the tangency 
requirement, a new approximation to the solution is obtained 
from Eq. (D27) , where the sign chosen for the radical is 
important, 

g. The iteration procedure continues with the new 

u (x, 0) and T until both quantities converge to a solution, 
w 

h. Once a solution is obtained, the frees tream Mach 
number for a given flow condition is calculated from 
Eq. (D18) and the local Mach nxamber distribution from 
Eq. (DIO) . 
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APPENDIX E 



HODOGRAPH TECHNIQUES FOR TRANSONIC FLOW 



The hodograph plane represents a coordinate transfor- 
mation in which the velocity components are the independent 
variables. Spatial coordinates x and y are replaced by 
velocity components ti and that is, the physical coordi- 
nates are represented by x(u,v) and y(u,v) [Ref. 26]. This 
coordinate transformation is convenient because it enables 
the simple presentation of data or solutions and, more 
importantly, it linearizes the two-dimensional planar tran- 
sonic flow equation that is nonlinear in the physical plane. 
This equation is known as the Tricomi or Tricomi-Euler Equa- 
tion. Boundary conditions are difficult to satisfy [Ref. 15] . 

The small perturbation transonic equation without heat 
addition can be written as: 




(El) 



where for irrotational flow 




0 



(E2) 



Using the hodograph transformation, Eqs. (El) and (E2) 



become : 
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(E3) 



and 




3v' 3u 



(E4) 



Eqs . (E3) and (E4) are the linear transonic hodograph 
equations . 

As is mentioned in the introduction, the linearized 
small perturbation equations for supersonic flow are hyper 
bolic and those for subsonic flow are elliptic. The 
transonic hodograph equations change from elliptic to 
hyperbolic when the coefficient of 3y'/9V changes sign; 



The critical velocity, where the change of sign occurs, 
can be directly solved from Eq. (E5) . 



Therefore, for u' < u'*, Eq. (E3) is elliptic, and when 
u' > u'*, it is hyperbolic. The features of transonic flow 
are captured by the transonic hodograph equations . Since 
the flow occuring in the external region and on the bounding 
streamtube is considered to be isentropic two-dimensional 



(1-M^) - M^(y+1)u' = 0 

00 CO 
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a.nd planar, the hodograph method gives a conceptual idea of 
''^hat the resulting flow configuration should be over an 
infinite wedge of prescribed shape. 

Guderley and Yoshihara [Ref. 15] use the flow angle 0 
and n, defined in Eq. (E7) , as the independent variables. 

n = (y+ 1) ( (w-w*) /w*) (E7) 

r~2 2 ^ 

where w = y/u + v . With this selection of independent 
variables, the streamf unction in the hodograph plane is; 

♦tin • *ee ' “ (^8) 

and the velocity potential function is; 

♦nn ■ " *ee = « 

The formulation is for a flow with constant stagnation 
temperature and, hence, can be applied to the solution for 
the bounding streamtube. For flow over a wedge, Guderley 
and Yoshihara found a steady solution for precisely unity. 

Fig. lEa shows the flow over the bounding streamtube 
in the physical plane. Points A, B, C, D, and E are shown 
mapped in the hodograph plane. Fig. lEb. A limiting Mach 
line is shown in the physical plane. Any characteristic in 
the region bounded by the sonic line and the limiting Mach 
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(b) Flow in the Kodograph Plane 



FIGURE IE. Sonic Flow over Bounding Streamtubes 
(Physical and Kodograph Plane 
Representations ) 
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line can interact with the sonic line. Any characteristic 
downstream of the limiting Mach line does not interact with 
the sonic line, and, hence, the flow is not affected by the 
subsonic flow upstream of the sonic line. 
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APPENDIX F 



DERIVATION OF THE SLUING RATES NECESSARY TO 
PRECLUDE SIGNIFICANT PHASE DISTORTIONS IN 
A LASER BEAM AT MACH 1.0 



When the phase difference between two points in a laser 
beam is greater than a quarter of a wavelength, the deviation 
of the beam from its original direction is sufficient to 
cause considerable distortion. The range of laser sluing 
rates necessary to prevent such phase distortions and, hence, 
to preclude defocusing of the beam at Mach 1.0 is derived. 

Consider a laser beam that is slued as shown in Fig. 1. 
The Mach nximber of the air relative to the beam varies along 
the beam as: 



M = 





(FI) 



All Mach nvimber regimes occur along the beam. If L repre- 
sents the radial distance between two stations on the beam 
where the Mach numbers are M^ and M 2 , it can be determined 
from Eq. (FI) that L is equal to: 



(M2~M^) a 

n 



(F2) 



The number of waves (N) of laser radiation in the region 
of length L is given by: 
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(F3) 



N = - = 

X c 

where X is the wavelength of the laser beam and is equal 
to nv/c. V (radians sec is the laser frequency equal 
to approximately 2 x 10^^ Hertz for a CO 2 laser, and c (m 
sec ) is the speed of light. Therefore, the change in the 
number of waves in the beam length L, due to a change in 
the index of refraction, is obtained by differentiating 
Eq. (F3) . 



dN = -VA^. 
c 



(F4) 



The change in the refractive index is related to the 
density perturbations in the beam, a known quantity, through 
the Gladstone-Dale Law [Ref. 37] , a simplified version of 
the Lorentz-Lorenz Law for gases with n = 1, Eq. (F5) : 



n = 1 + kp , (F5) 

where k is the Gladstone-Dale constant which is approximately 
equal to 2.25 x lo”^ m^kgm for air with wavelengths beyond 
the near infrared (X > ly) • Differentiating Eq. (F5) and 
substituting it and Eq. (F2) into Eq. (F4) , the desired 
expression for the change in the number of waves is obtained 
as : 



DN 



k(M--M, ) avp^ 

( ) 

fic 






(F6) 
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It is of interest to determine the range of sluing 
rates for which the change in the number of waves is less 
than a quarter of a wavelength at Mach 1.0 from Eq. (F7) ; 



4k (M-,-M. ) avp 

> L-1 



00 




(F7) 



c 




and is the minimum (maximum negative) density perturba- 

tion in the beam at Mach 1.0/ points A and B in Fig. 11, 
respectively. As an example, using M 2 = 1.001, M^ = 0.999, 
the flow properties in Table I, and the maximum difference 

in the density perturbations at Mach 1.0 from Fig. 11 of 

-4 -3 

1.5 X 10 , the sluing rate must be greater than 7.5 x 10 

radians/sec . 
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APPENDIX G 



SUBSONIC AND SUPERSONIC THERMAL BLOOMING 
AND COMPUTER PROGRAM (BLOOM) 

Tsien and Beilock [Ref. 38] have obtained solutions for 
the linearized equations of motion with heat addition for 
subsonic and supersonic flow. The solutions are for a line 
heat source q in a uniform, two-dimensional planar flow 
which can be considered to be a perturbation on an initially 
one-dimensional flow. 

In 1971 and 1972 studies were conducted at NPS by Fuhs 
[Ref. 39] on the applicatin of external heat to the reduction 
of drag. These studies required thorough analysis of heat 
addition in both subsonic and supersonic flows. The analysis 
of the flow with heat addition found two other applications. 
One application. was external burning assisted projectile. This 
application was studied by CDR. W. J. H. Smithey [Ref. 40] . 

The other application was recognized by Fuhs [Ref. 4]; this 
involves the density inhomogeneity arising from quantum 
inefficiency. Professors Biblarz and Fuhs [Ref. 5] developed 
the capability and extended the analysis to predict density 
inhomogeneity in a large gas dynamic laser; see the paper 
by Fuhs, Biblarz, Cawthra and Campbell [Ref. 6] for details 
of the experiment. Subsequent to the application of linear- 
ized solutions to flow internal to a laser, it became apparent 
that the analysis could predict flow properties in thermal 
blooming. In a series of papers and presentations, Fuhs, 
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Biblarz, Burden and Carey [Refs. 1 , 8 and 41-43] developed 
the appropriate analysis for linearized thermal blooming. 

1. SUBSONIC FLOW WITH HEAT ADDITION 

The perturbations to the flow [Ref. 38] are given by: 



a* = 5L 



V = 



p' = - 



p' = - 



(x2 + eV) 

(Y-l) 6q y 

2 ^YP„u„ V 2 ^ o 2 2 , 

^ 0 ° “ (x + 3 y ) 

(Y-l)M,„q 



2 2 2^ 
oo-iToo + By) 



2TTa„3p 
(Y-l)M„q ^ 

(x^ + bV) 



(Gl) 



(G2) 



(G3) 



- 6(y) i(x) (G4) 



3 

Map 
00^00 ^00 



where the symbols are defined in the list of symbols and 

/ ^ 

3 = /I - M for subsonic Mach numbers. The second term in 
Eg. (G4) is the wake of the line heat source q. The 
temperature perturbation can be obtained from the equation 
of state. 

Biblarz and Fuhs [Ref. 6] in applying Eqs. (Gl) through 
(G4) to laser internal aerodynamics have presented the 
integral equations used to solve the unbounded heat release 
region problem (i.e. subsonic thermal blooming) and are 
listed below for completeness. Since the heat source is 
volume distributed, instead of a single infinitesimal rod, 
an integral representation was used to account for all 
contributions to the flow properties over the volume. The 
geometry involved in the calculation is shown in Fig. IG. 
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FIGURE IG. Geometry of Flow 



v' 



p' = - 



_ (y-1) 


/ / 
5 n 


(x-C) hCC/T,) 


dCdn 


2"BP„u„ 


(x-c)^ + e^(y-n) 


_ (y-1)B 


/ / 
C n 


(y-n) h(c,n) 


d^dn 


2"YP„p„ 


(x-c)2 + 8^(y-p)^ 


(y-1)m. 


L / / 


(X-;) h(s,7i) 





5 n ,^.5)2 ^ 



p' 



27ra^^Sp^5 n (x-?)^ + 6^(y-n) 



/ / . dcdn 



I j h(§,n) 6(y-n) I(x-0 d^dn 



a u p f n 
00 00^00 ^ ' 



(G5) 

(G6) 

(G7) 

(G8) 
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The heat quantity h(x,y) is related to the laser output 
intensity through the kinetic lag that may be present. To 
develop an equation incorporating the lag between laser 
energy and heating, it is necessary to examine the kinetics 
of relaxation. An exponential function describes the 
relaxation with a relaxation time x [Refs. 44 and 45]. The 
heat release function is given by: 



h(x,y) = I^a / f(x',y) exp[- ] 



dx’ 



(G9) 



where Q = I^a f(x,y) is the laser energy per time per 
volume. The volume can be interpreted as an area times the 
span of the heat addition zone across the planar flow. 

If significant variations of beam intensity in the radial 
(beam) direction exist, the two-dimensional heat addition 
function can be approximated from the three-dimensional 
h(x,y,r) by an average: 



^i+1 



h(x,y) 



= J. / 



Ar 



h(x,y,r) dr 



(GIO) 



r . 

1 



where Ar = 



The beam is divided into segments Ar 
in length and Eqs. (G5) through (G8) are solved for each 
segment. 

If the heat release h(x,y) is constant, the equations 
can be integrated in closed form. The resulting integrals 
are: 
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(Y-Dl^a 




27TYP„3u^ M(x,y;g,n,3) 


(Gll) 


(Y-1) I a 




'' ■ 2typ,6u„ N(x,y;5,n,6) 


(G12) 


(Y-l)MI a 




P ~ 2i'a 6p M(x,y;g,n,3) 

00 ^00 


(G13) 



(Y-l)M^I a 

P = - 3— ^ M(x,y;C,n,6) 

(Y-1) l a 

“ — 2 // <5(7-11) I(x-C) dCdn (G14) 

^00 ^00 Poo 5 ri 



where 

M(x,y; C/n, 3 ) 



(y-n) 

2 



In 



(x-g)^ + g^(y-n)^ 
6^(y-n)^ 



+ 



(x-g) 

B 



tan 



-1 



3 (y-n) 

(x-g)J 


■ ^2 


- 


■ "^1. 



^2 

«1 



N(x,y; g,n,6) = 



(x-g) 

2 

26 ^ 



In 



(x-g)^ + 62(y_^)2 
3 ^ 



- tan-1 

^2 6 



(x-g) 
3 (y-n) 



ho ' 



The first term in Eqs. (G8) and (G 14 ) correspond to the 
density perturbations in regions outside the wake and the 
last term represents the contribution of the wake and can 
be interpreted as the consequence of the heat convected 
along the streamlines by the flow. 
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2. SUPERSONIC FLOW WITH HEAT ADDITION 

The perturbations to the flow are also formulated by 
Tsein and Beilock [Ref. 38] as: 



u' = - 6(x - 3 y) 



(G15) 






(G16) 



(Y-l)M^q 



(G17) 



p' = - 



(Y-l)M^q 

2a ^3p 
*•^00 '^^00 



6(x - 3y) - — 6 (y) l(x) 



M a P 
^*00 ”*00 ^00 



(G18) 



Fuhs [Refs. 4 and 39] has modified these equations into 
an integral form as: 






(G19) 



' r'oo 00 Q 

(Y-I)M^ s 

P' = 2a"Tb ~~ ^ h(C,n) sin U ds 

^ “oo ^ ^00 



(G20) 



(G21) 



(Y-I)M^ S 

p' = ^ / h(5,n) sin y dS 



2a„ 3 p„ 0 



/ / h(?,n) 6(y-n) i(x-5) d5dn (G22) 



M„a„ p„ c n 
00 00 00 



where the synibology is the sane as that used in the subsonic 
equations except that B = - 1 for supersonic Mach numbers, 
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3. COMPUTER PROGRAM 



The FORTRAN computer program BLOOM, included at the end 
of this appendix, was written to compute and/or plot the flow 
perturbations (u', v', p' and p') in the subsonic and super- 
sonic flow with heat addition. Using a modified version of the 
computer program DIDER [Ref. 5] as a subroutine, the subsonic 
perturbation equations, Eqs . (G5) through (G9) , and the 

supersonic perturbation equations, Eqs. (G19) through (G22) , 
are numerically integrated to obtain the linear thermal 
blooming solutions. The program calculates the two- 
dimensional field at any cross section of the laser beam 
where the energy release distribution is known. Finite 
kinetics are incorporated as described by Eq. (G9) . 

The program (BLOOM) is versatile in that it can calculate 
the flow quantities for several subsonic and/or supersonic 
freestream Mach numbers, heat release distributions, relaxation 
rates, and freestream flow properties in a single computer run. 
Since the length of time for each individual case depends 
primarily on the size of the flow region and on the mesh that 
is superimposed upon it, time requirements must be determined 
on an individual basis. For the results presented in Figs. 4 
through 7 and Figs. 12 and 13 from the program BLOOM, the 
two-dim.ensional region of interest was 60 cm. square centered 
at the origin of the Gaussian heat release distribution and 
was descretized into a mesh of cells 1 cm. square for 
subsonic flow (M = 0.999) and into rectangular cells whose 

oo 

size is determined by the Mach angle for supersonic flow 
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(^00 “ 1.001). Supersonically, the ratio of cell height Ay 
to cell length Ax is equal to tan y. The subsonic results 
take approximately 25 minutes and 140 K of memory (without 
plotting subroutine CONTUR) while the supersonic results 
take approximately 3 minutes and 285 K of memory (without 
plotting routine CONTUR) to compile and execute using an 
IBM 360 computer. 

The program can also calculate constant heat release 
distributions for single and multi-pass laser beam configura- 
tions with or without wall reflections and kinetics. In its 
present form, the program can not calculate supercritical 
(mixed) flows across the beam. 

A main feature of the supersonic portion of the program 
is that it sets up the characteristic lines belonging to each 
individual cell (i.e. their diagonals) and adds up the 
contributions for each flow property from the characteristics 
and, in the case of the density perturbation, from the wake 
at each cell. The calculations along the characteristics 
are facilitated by the choice of cell shape and size as 
described above. 
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CALCULATF CONTRIBUTIONS FROM CHARACTERISTICS FOR SUPERSONIC FLOW 
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